-
-
Notifications
You must be signed in to change notification settings - Fork 5.5k
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Improve the Type consistency of Special Functions (Fixes #17474) #18584
Conversation
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
note the whitespace error
base/special/gamma.jl
Outdated
x = float(z) | ||
t = Int(s) | ||
if typeof(x) === typeof(z) && typeof(t) === typeof(s) | ||
throw MethodError(zeta,(s,t)) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
throw(...)
and why is this a MethodError exactly? could use a comment
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Thanks.
It is a method error because this is a fallback definition.
If the most specific implementation of this is foo(::Number)
,
and if once you convert it to a float
, the most specific implementation is still foo(::Number)
,
because the type hasn't changed,
then we don't have a working implementation for this type.
This is the case for trigamma(::BigFloat)
I'll try and put that into a comment.
Should I write a block comment before all the fallback methods?
or write one per throw
?
base/special/gamma.jl
Outdated
function invdigamma(y::Float64) | ||
# Closed form initial estimates | ||
# Implementation of fixed point algorithm described in | ||
# "Estimating a Dirichlet distribution" by Thomas P. Minka, 2000 |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
don't change the indent of comments
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I moved the comment into the inside of the function it was commenting; So it would not obscure the doc-string.
Should I not indent comments, that are inside functions? (honest question, I'm not sure on the convention)
I kinda think maybe I should have just added that comment to the doc-string.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
if it's inside a function then it should be indented to the same level that it would be if it were code on the same line, generally - so a multiple of 4 spaces here
test/math.jl
Outdated
@test typeof(digamma(big"2")) == BigFloat | ||
@test_throws MethodError trigamma(big"2") | ||
@test_throws MethodError trigamma(big"2.0") | ||
@test_throws MethodError invdigamma(big"2") |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
are these not available in mpfr/gmp?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Correct. They are not available, or if they are, then we are not importing them.
Let me check mpfr docs
Edit:
They are not available in MPFR:
http://www.mpfr.org/mpfr-current/mpfr.html#index-Special-functions
GMP doesn't have it since (as I understand it) these functions almost never return integers.
I need to make test cases for to make sure they throw method errors, or work correctly. And for my future reference: the test cases can use: julia> eta(Complex(big"52.0")) |> real |> BigFloat #A
9.999999999999997779553950749686919152736663818359375000000000000000000000000000e-01
julia> eta(big"52.0") #B
9.999999999999997779553952297414836228661208246910286674974797896945379409710455e-01
julia> eta(Complex(52.0)) |> real |> BigFloat #C
9.999999999999997779553950749686919152736663818359375000000000000000000000000000e-01 It should be A==B not C, or a Method Error Similarly: julia> zeta(Complex(big"52.0")) |> real #A
1.000000000000000222044604925031308084726333618164062500000000000000000000000000
julia> zeta((big"52.0")) #B
1.000000000000000222044605079804198399932009420465396423665432943893439236654051
julia> zeta(Complex(52.0)) |> real |> BigFloat #C
1.000000000000000222044604925031308084726333618164062500000000000000000000000000 Again A==B not C, or a Method Error |
Thinking exactly how I want everything to behave. The functions in
|
base/special/gamma.jl
Outdated
# "Estimating a Dirichlet distribution" by Thomas P. Minka, 2000 | ||
for f in (:zeta, :eta) | ||
@eval begin | ||
$f{T<:Union{Base.BitInteger,Float32,Float16}}(x::t) = oftype(float(x), $f(Complex128(z))) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Complex128(x) ?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Complex128
is a type alias for Complex{Float64}
.
I guess cos it is made of 2 Float64s
It is odd that it is converting a real to a complex though.
Looking at it, it is because that line doesn't even compile. (I got distracted when writing this function i guess).
I know I've fixed it already in my local copy.
x::t
is x::Complex{T}
Or probably x::Complex{Union{Base.BitInteger,Float32,Float16}}
Can't remember
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
sorry, I was very cryptic, I meant z should be changed to x in Complex128(z) :) but of course you already fixed that
Sorting this was a bit more complex that I thought. I have moved it from WIP into RFC. |
base/special/gamma.jl
Outdated
ofpromotedtype(as::Tuple, c) = convert(promote_type(typeof.(as)...), c) | ||
ComplexOrRealUnion(TS...) = Union{(ComplexOrReal{T} for T in TS)...} | ||
|
||
const types_le_Float64 = (Float16, Float32, Float64, Base.BitInteger.types...) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Is this a good name for this?
It is the collection of types that are less-than or equal to precise than Float64, at representing Float64 inputs.
I don't think this follows naming convention
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I'd probably just use the tuple in the loops and avoid defining the variable. It's only used three times.
base/special/gamma.jl
Outdated
# they also implement there own versions of the functions | ||
|
||
|
||
ofpromotedtype(as::Tuple, c) = convert(promote_type(typeof.(as)...), c) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Should this move to promotion.jl
?
It is the merging of oftype
with promote
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
If you use the strategy I outlined (which might not work) then I don't think this function is necessary. I doubt that it will find much use outside this function anyway, so I think it can stay here if it ends up being useful.
base/special/gamma.jl
Outdated
|
||
|
||
ofpromotedtype(as::Tuple, c) = convert(promote_type(typeof.(as)...), c) | ||
ComplexOrRealUnion(TS...) = Union{(ComplexOrReal{T} for T in TS)...} |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I would rather this was ComplexOrRealUnion{TS...}
or something like that.
But I can't seem to get typealias
to take a variable number of arguements
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Maybe split the definitions into two line to avoid this function. It might be easier to read.
if typeof(t) === typeof(s) && typeof(x) === typeof(z) | ||
# There is nothing to fallback to, since this didn't work | ||
throw(MethodError(zeta,(s,z))) | ||
end |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I think this should have a second fallback to zeta(::Number, Number)
maybe.
Or possibly to prevent looping just try and convert s
to a float and see how that goes.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
See my other comment about the zeta
.
Is there anything I need to do now, for this to proceed? |
base/special/gamma.jl
Outdated
@@ -106,7 +108,8 @@ function lgamma(z::Complex{Float64}) | |||
-2.2315475845357937976132853e-04,9.9457512781808533714662972e-05, | |||
-4.4926236738133141700224489e-05,2.0507212775670691553131246e-05) | |||
end | |||
# use recurrence relation lgamma(z) = lgamma(z+1) - log(z) to shift to x > 7 for asymptotic series | |||
# use recurrence relation lgamma(z) = lgamma(z+1) - log(z ) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
unnecessary space in log(z )
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
fixed, thanks
@tkelman Is there anything I need to do, that is holding this back? |
someone who knows this code should respond to your inline questions |
@andreasnoack you wrote a lot of this 2 years back; right? |
base/special/gamma.jl
Outdated
# and results converted back. Similar for BitIntegers (eg Int32), but no converting back | ||
# Otherwise, we need to make things use their own `float` converting methods | ||
# and in those cases, we do not convert back either as we assume | ||
# they also implement there own versions of the functions |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
their own
base/special/gamma.jl
Outdated
# Otherwise, we need to make things use their own `float` converting methods | ||
# and in those cases, we do not convert back either as we assume | ||
# they also implement their own versions of the functions, with the correct return types. | ||
# Otherwise, if they do not implement their version of the functions we |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
whitespace error
whitespace was fixed, implementation still needs review though
@simonbyrne @stevengj this basically implements what we discussed in #17474 |
test/math.jl
Outdated
|
||
@test Base.Math.f64(complex(1f0,1f0)) == complex(1.0, 1.0) | ||
@test Base.Math.f64(1f0) == 1.0 | ||
@test Base.Math.ofpromotedtype((complex(1f0, 1f0), 1.0), 5.0) == complex(5.0) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
should these check types too with ===
?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Yes, they should. (I will change)
base/special/gamma.jl
Outdated
oftype(z, polygamma(m, f64(z))) | ||
end | ||
|
||
for T1 in types_le_Float64, T2 in types_le_Float64 |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I think meta programming sometimes can be a bit convoluted to read. Would it be simpler to define something like
for T in (Float16,Float32,Complex{Float16},Complex{Float32})
zeta(x::T, y::T) = T(zeta(f64(x),f64(y)))`
end
zeta(x::Integer,y::Integer) = zeta(f64(x),f64(y))
zeta(x::Number, y::Number) = zeta(promote(x,y)...)
zeta(x::BigFloat,y::BigFloat) = error
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This is indeed a particularly convoluted piece of metaprogramming.
I rewrote it about 3 times in different ways.
but I think I'll try again.
Get rid of some of the if T1<:Integer && T2<:Integer
by using a separate loop for ints.
I have been trying to avoid things like
zeta(x::BigFloat,y::BigFloat) = error...
, because that will mean that if someone does write a zeta(::BigFloat, ::BigFloat)
it will have issues with overwriting a function name. And will trigger #265 problems.
In general, rather than exclude the types in Base
that a function does not work for,
I am trying to include only the types in Base which a function does work for.
Which should make extending with custom number types from packages sane.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Since float(::BitInt)::BigFloat
, wouldn't the version below accomplish what you want?
for T in (Float16,Float32,Complex{Float16},Complex{Float32})
zeta(x::T, y::T) = T(zeta(f64(x),f64(y)))`
end
zeta(x::Integer,y::Integer) = zeta(float(x),float(y))
zeta(x::Number, y::Number) = zeta(promote(x,y)...)
zeta(x::BigFloat,y::BigFloat) = error
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Didn't see your response before I made mine here
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I think the problem is that a line like zeta(x::Number, y::Number) = zeta(promote(x,y)...)
is really simple and easy to read relative to the implementation suggested in the PR. The problem with the version I suggest is, of course, that it is an infinite recursion if no promotion happens. Therefore the explicit error for BigFloat
. Rumors are that #265 is close to getting fixed so maybe the redefinition concern shouldn't weight too much.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Did you see my implementation of: function zeta(s::Number, z::Number)
etc?
Line 694.
This doesn't suffer from infinite recursion.
I might be able to delete quiet a few definitions, and just fall back to it.
(It might need a promote
added after the float
)
It might end up a lot like yours, I think.
But avoid the explict erroring for BigFloats.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
That could work. If you restrict the method in 694 to promote to the same element type for the two arguments. Then the method here can be restricted to a single loop instead of looping through the combinations of both argument types.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Can zeta's first argument just always be a Float64?
I think it can.
I don't think there is a meaningful difference to giving it Int
(maybe a microspeed up).
If it can't then promoting to the same type won't work.
But I think it can,
so that is find, and Integers can be treated as just numbers, when in the first arg position (as well as in the second).
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Yes. I think it is fine to remove Int
from the signature.
base/special/gamma.jl
Outdated
# they also implement there own versions of the functions | ||
|
||
|
||
ofpromotedtype(as::Tuple, c) = convert(promote_type(typeof.(as)...), c) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
If you use the strategy I outlined (which might not work) then I don't think this function is necessary. I doubt that it will find much use outside this function anyway, so I think it can stay here if it ends up being useful.
base/special/gamma.jl
Outdated
|
||
|
||
ofpromotedtype(as::Tuple, c) = convert(promote_type(typeof.(as)...), c) | ||
ComplexOrRealUnion(TS...) = Union{(ComplexOrReal{T} for T in TS)...} |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Maybe split the definitions into two line to avoid this function. It might be easier to read.
base/special/gamma.jl
Outdated
ofpromotedtype(as::Tuple, c) = convert(promote_type(typeof.(as)...), c) | ||
ComplexOrRealUnion(TS...) = Union{(ComplexOrReal{T} for T in TS)...} | ||
|
||
const types_le_Float64 = (Float16, Float32, Float64, Base.BitInteger.types...) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I'd probably just use the tuple in the loops and avoid defining the variable. It's only used three times.
base/special/gamma.jl
Outdated
|
||
for f in (:digamma, :trigamma, :zeta, :eta, :invdigamma) | ||
@eval begin | ||
$f(z::ComplexOrRealUnion(Base.BitInteger.types...)) = $f(f64(z)) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Are we sure that it is so dangerous to convert all integers to Float64
?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
The problem is not accurasy (for real Ints at least), it is inconsistant return types.
Loss of accuracy does not occur until beyond ints 2^53 -- which I'm yet to find a case that it matters for any of these functions, as they tend to hit an asymptote long before then
We want BitInteger
(Int64
etc) to have Float64
results from these functions (at least that is current behavior).
The main deal is we want BigInt
s to have BigFloat
results if possible,
and to throw method errors if not (rather than having Float64
results).
Which is most of #17474
Eg right now
julia> zeta(big"52")
1.0000000000000002
julia> zeta(big"52.0")
1.000000000000000222044605079804198399932009420465396423665432943893439236654051
We really want to avoid reporting BigFloat
results, that were calculated with Float64
's behind the scenes, and then converted. as in the examples in #18584 (comment)
I don't think that kinda thing ever happens for BigInt
, right now.
Because right now most BigIntegers
functions return Float64s
Even if there is more accurate BigFloat
method available in MPFR
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Got it. Thanks for explaining. See my updated "proposal" below.
if typeof(t) === typeof(s) && typeof(x) === typeof(z) | ||
# There is nothing to fallback to, since this didn't work | ||
throw(MethodError(zeta,(s,z))) | ||
end |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
See my other comment about the zeta
.
@oxinabox Sorry for submitting my comments so late. I made them a while ago and thought they were visible. Still adapting to the review system. |
Right, sorry for the delay. I think I have made most of @andreasnoack 's suggestions that we ended up with.
Also now the return type of |
base/special/gamma.jl
Outdated
# It is safe to convert `s` to a float as underflow will occur before precision issues | ||
zeta(s::Integer, z::Number) = zeta(oftype(z, s), z) | ||
function zeta(s::Integer, z::Integer) | ||
x=float(z) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
tabs -> spaces
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
ah, bother. I redid my vimrc, and must have forgotten to change that back.
Thanks
base/special/gamma.jl
Outdated
ζ = zeta(f64(s), f64(z)) | ||
|
||
if typeof(s) <: Complex || typeof(z) <: Complex | ||
convert(Complex{$T2}, ζ) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Shouldn't it be promote_type($T1,$T2}
instead of of $T2
?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
It used to be, but after thinking a bit more about zeta
,
I started to think that it should maybe be driven only by the type $T2
.
But I now think I was going down the wrong line.
- I wanted
zeta(::BigInteger, ::Float64) :: Float64
(which is what depending only on$T2
gives me) promote_type($T1,$T2}
, gives mezeta(::BigInteger, ::Float64) :: BigFloat
-- which is a lie as that precision is not there.- I now think that it should just
MethodError
, which I can get by deleting a definition (or two)
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
+1 to MethodError
if either argument is Big
.
|
I think this applies to the current version as well. There you'll get a To summarize, I think there are three solutions
I recognize the problem with |
Note that we do manually Line 452 in 68af05e
It's not pretty, but I think it is the least worst option. |
I wonder if |
#17057 will include that flag :) |
base/special/gamma.jl
Outdated
Compute the polygamma function of order `m` of argument `x` (the `(m+1)th` derivative of the | ||
logarithm of `gamma(x)`) | ||
Compute the polygamma function of order `m` of argument `x` | ||
(the `(m+1)th` derivative of the logarithm of `gamma(x)`) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
"th" shouldn't be code highlighted
the other trick is to have a different function, e.g.
|
bump |
I think this is still good to go. Unless the suggestion from @simonbyrne in #18584 Rebasing is getting harder the more time passes. |
Thanks for the reminders. Especially relevant as there are few open PR's that move this code around but aren't intending to change behavior. cc @ararslan, what do we want to do here short-term? |
@tkelman Tough call, as obviously a lot of work has gone into this PR and I've also put a lot of work into moving these things out. I guess it depends on whether getting this stuff out for 0.6 is a priority. If so, I'd say this should be reopened against SpecialFunctions (a nontrivial task), otherwise I'd say merge this and I can update things accordingly. |
Since the package is now registered and #20427 is just about ready, I'm going to assign @ararslan as a next-step to port this over to the package version. Thanks for the contribution and sorry this didn't get merged a few months ago, we'll get these improvements into the package (where people even on release Julia can use it right away). |
As I found in #17474,
the various special functions tend to just blindly cast everything to Float64.
Regardless of if it makes sense.
This is my attempt to fix that.
And while I am at it, to make the return types also a little more consistent across the functions.
This is WIP, I'm not actually sure my tests pass, but I would like feedback as to if I am going about this in a sane way.Though these do need to move to SpecialFunctions.jl
(cf: JuliaMath/SpecialFunctions.jl#12)